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ABSTRACT 

We study cell count moments up to fifth order of the distributions of haloes, of 
halo substructures as a proxy for galaxies, and of mass in the context of the halo model 
and compare theoretical predictions to the results of numerical simulations. On scales 
larger than the size of the largest cluster, we present a simple point cluster model in 
which results depend only on cluster-cluster correlations and on the distribution of the 
number of objects within a cluster, or cluster occupancy. The point cluster model leads 
to expressions for moments of galaxy counts in which the volume-averaged moments 
on large scales approach those of the halo distribution and on smaller scales exhibit 
hierarchical clustering with amplitudes Sk determined by moments of the occupancy 
distribution. In this limit, the halo model predictions are purely combinatoric, and have 
no dependence on halo profile, concentration parameter, or potential asphericity. The 
full halo model introduces only two additional effects: on large scales, haloes of different 
mass have different clustering strengths, introducing relative bias parameters; and on 
the smallest scales, halo structure is resolved and details of the halo profile become 
important, introducing shape-dependent form factors. Because of differences between 
discrete and continuous statistics, the hierarchical amplitudes for galaxies and for 
mass behave differently on small scales even if galaxy number is exactly proportional 
to mass, a difference that is not necessarily well described in terms of bias. 
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1 INTRODUCTION 

Describing the properties of the distribution of matter in 
the universe in terms of the masses, spatial distribution, 
and shapes of clusters, or haloes, is an enterprise with a 
long history (Neyman & Scott 1952; McClelland & Silk 
1977; Peebles 1980; Scherrer & Bertschinger 1991; Sheth & 
Saslaw 1994; Sheth 1996b). Recently, with the new ingredi- 
ent of a universal halo profile found in numerical simulations 
(Navarro, Frenk, & White 1997; Moore et al. 1999; Navarro 
et al. 2004), interest in the model has been rekindled (Seljak 
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2000; Ma & Fry 2000b; Peacock & Smith 2000; Scoccimarro, 
Sheth, Hui, & Jain 2001). This model is not seen as literally 
true, but its constructions give plausible estimates for cor- 
relation functions because at a given scale, density-weighted 
statistics are dominated by the highest density systems, the 
collapsed haloes. The model has been shown to reproduce 
two-point and higher order density correlation functions in 
simulations, and, with a carefully chosen halo mass func- 
tion and 'concentration parameter,' can be consistent with 
self-similar stable clustering (Ma & Fry 2000a; Smith et al. 
2003). Among its many other applications to weak gravi- 
tational lensing, pair velocities, the Ly-a forest, and CMB 
foregrounds, we find that the halo model also allows us to 
address the different behaviors of the continuous mass den- 
sity and of discrete objects such as galaxies. 
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In this paper we reexamine the statistical behavior of 
integral moments of total mass or number counts in the 
context of the halo model. Our results, formulated directly 
in position space, complement and extend those of Scocci- 
marro, Sheth, Hui, & Jain (2001) in the Fourier domain. In 
Section 2 we present definitions of the various statistics we 
need and introduce generating function tools that will be 
applied in later sections. In Section 3 we apply the probabil- 
ity generating function machinery for a system of identical 
clusters in the point cluster limit, a model we call the 'naive 
halo model,' to express the statistics of counts in cells in 
terms of properties of the halo number and halo occupancy 
distributions. In Section 4 we compare the model to results 
obtained from numerical simulations. We find that the naive 
point cluster model describes the qualitative behavior but 
fails in quantitative detail, but insight gathered from the 
model in the generating function formalism is easily applied 
in the full halo model. This leads us in Section 5 to consider 
the halo model in its full detail, summing over haloes of 
different mass, with both halo occupations and halo correla- 
tions functions of halo mass. Properly interpreted, the naive 
point cluster results obtained using the generating function 
continue to apply when averaged over halo masses and over 
galaxy positions within a halo. This allows us to extend to 
small scales, where haloes are resolved, introducing geomet- 
ric form factors for haloes that can no longer be consid- 
ered as points. Working directly in space instead of in the 
Fourier transform domain allows us to exhibit manifestly 
symmetries under particle exchange at all orders; avoids the 
necessity to introduce an approximate factoring of window 
function products W1W2W12 ~ WiW£, etc.; and avoids the 
necessity to make any assumptions or approximations about 
configuration dependence. Known forms of the halo mass 
function, bias, and occupation number allow us to compute 
from first principles results in scale-free and specific cos- 
mological models. Section 6 contains a final discussion, and 
some technical results are included in appendices. 



2 STATISTICAL DEFINITIONS 

We study statistics of the continuous mass density and of 
discrete objects that for convenience we denote as "galax- 
ies." For a galaxy number distribution that is a random sam- 
pling of a process with a smooth underlying number density 
field n(r), factorial moments of the number of galaxies in a 
randomly placed volume directly reflect moments of the un- 
derlying continuous density field (Szapudi & Szalay 1993), 

(w [fcI ) = #*#*, (1) 

where N [k] = N\/(N - k)\ = N(N - 1) • • ■ (N - k + 1), 
and the moments p. k are volume averages of corresponding 
moments of the underlying density field, 



^k = yj: d 3 ri---d 3 r fc ^ fc (ri,...,r fc ), 
(n(r-i) •■■n(r k )) = n k n k (n, . . . ,r k ), 



(2) 
(3) 



typically integrated over a spherical volume of radius R. 
Moments of powers ( N ) then contain contributions arising 
from discreteness; for k = 2 through 5 these are 



(N 3 )=N + 3N 2 ft 2 +N 3 fi 3 , (5) 

(iV 4 ) =7V + 7iV 2 /i 2 + 6iV 3 ^3 + iV 4 /i4, (6) 

(N 5 )=N + 15iV 2 /l2 + 25N 3 fj, 3 + WN 4 fJ,4. + N 5 fi 5 . (7) 

In the limit N — ( N \ ^> 1 the highest power of N domi- 
nates and ( N k \ = N k jl k , as for a continuous density; the 
factorial moment, in removing the lower order or discrete- 
ness terms, leaves a discreteness corrected moment that re- 
flects only spatial clustering. The moments p, k can be addi- 
tionally separated into irreducible contributions 6, as 



M2 = 1+6, 

#3 = 1+362+13, 

#4 = 1 + 6C2 + 3|| + +46 + 6, 



(8) 

(9) 

(10) 



#5 = 1 + 10C2 + 106 + 156 + 1066 + 56 + 6, (11) 

also written as "connected" moments, 



(N [k] ) c = N k tk- 



(12) 



The relations between the fi k and the 6 can be summarized 
in the generating functions M(t) and K(t) = log M(t) (Fry 
1985), 






(13) 



M(t) = J2±N k !, k t k , 

k =o ' k=l 

With the factors 1/fc!, M and K are sometimes called ex- 
ponential generating functions of the moments JX k , 6- It is 
often found that the correlations vary with scale roughly as 
6 oc £ fe_1 . The hierarchical amplitudes S k are defined by 



6 — Sk £ 



(14) 



The normalization of 6 to S k then removes much of the 
dependence of 6 on scale. 

Generating functions provide an interesting connection 
between discrete and continuous processes. For a continuous 
variate x with moment generating function M c (t) = ( e tx >, 
the generator of a distribution of discrete counts N for which 
x is the local density is M d (t) = M c (e* - 1) (Fry 1985). This 
relation of generating functions provides directly the dis- 
creteness terms in equations (4-7). For discrete counts with 
probabilities Pjv, also useful is the probability generating 
function 



G(z) = J2 p ' 



N Z 



(15) 



( N 2 ) = N + N' 



1 12, 



(4) 



For a discrete realization of an underlying continuous num- 
ber density, G(z) is related to the exponential generating 
function of factorial moments by M(t) = G(t + 1) (Fry 1985; 
Szapudi & Szalay 1993). 



3 CELL COUNTS ON LARGE SCALES: THE 
POINT CLUSTER MODEL 

Using the tools introduced in the previous section we can 
now construct the generating function of total number count 
in the point cluster limit. On large scales, we expect that 
we can consider relatively compact clusters in their entirety 
to be either inside or outside of V. The total number of 
galaxies in a volume is then the sum over all the clusters in 
the volume, 
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N h 

* = £>, 



(16) 



where the number of clusters Nh and the number of galaxies 
Ni in each cluster are chosen randomly and at first we take 
the cluster occupation numbers Ni to be independent and 
identically distributed. A similar sum over clusters arises in 
situations ranging from the distribution of particle multiplic- 
ities in hadron collisions at high energy accelerators (Finkcl- 
stcin 1988; Hegyi 1994; Tchikilev 1999) to the distribution 
of rainfall totals (Rodriguez-Iturbe et al. 1987; Cowpertwait 
1994; Evin & Favre 2008). 

We can characterize the net count distribution directly 
for small counts and in general using the generating function 
G(z). Let p n be the probability of V containing n clusters, 
and let q„ be the probability that a cluster has n members. 
Because a cluster with no members is uninteresting, for sim- 
plicity we take qo = here; for the case qo 7^ 0, see Appendix 
A. Then, to have no count, we must have no clusters; to have 
total count 1, we must have one cluster with one member; to 
have total count 2 we can have one cluster with 2 members 
or two clusters with 1 member, and so on. The first several 
probabilities Pn that V contains N total galaxies are then 



Po = po 

Pi =Pi<7i 

Pi =Piq2 +P2ql 

P3 = PK73 + 2p 2 qiq2 + P?,qt 



(17) 
(18) 
(19) 
(20) 



P 4 =Pl<74 +P2(2<7l(j3 +<?2) +3p3<3 , l<3 , 2 +P4<7l (21) 

Ps = Pi <75 + P2(2<7iq 4 + 2g 2 Q3) 

(22) 



+p 3 (3<fr<73 + 'Sqiqi) + 4p4<7?g 2 + psqt- 



For these first few low order terms, Pn is the coefficient 
of z N in the composition of generating functions g^ [gi(z)\. 
This is the general result, as can be seen easily from the 
generating function for the total count probabilities Pn, 

G(Z) = <( z ^+N 2+ -+N Nh ^ = ^ z N t y h j 

= (M*)] Nh )=9HM*)]- (23) 

where the double angle brackets indicate an average over 
both the distributions of cluster occupancy and cluster num- 
ber (see Szapudi & Szalay 1993). 

We can also compute moments directly and using gen- 
crating functions. The mean of N — y] Ni is the number of 
haloes times the average occupation per halo, 



(N) = (N h Ni)=N h Ni. 



(24) 



The square N 2 — yj Ni yj Nj contains Nh terms with i — j 
and Nh(Nj — 1) terms with i 7^ j, 

( A 2 ) = ( N h N 2 + N h (N h - l)NiNj ) 

= Nh(Ni+N?ih,i) + NJl(l + $2, h )N? (25) 

Similar direct calculations give 

( AT 3 ) = ( N h N? + N h (N h - 1) 3N?Nj 

+ N h (N h - l)(N h - 2)N i N j N k ) (26) 

( A 4 ) = ( A h iV 4 + N h (N h - 1) (4iVf A, + 3A 2 N 2 ) 

+ N h (N h - l)(N h - 2) &N 2 NjN k 

+ N h (N h - l)(N h - 2)(N h - 3) NiNjNM ) (27) 



( A 5 ) = ( N h N? + N h (N h - 1) {5N?Nj + 10N?N 2 ) 

+ N l h 3] (WNfNjNk + 15N 2 N 2 N k ) 

+ N [ * ] WNfN-jNM + N l h 5] NiNjNMNrn ), (28) 

and from these, the discreteness corrected, connected mo- 
ments of total count are 

fi.2,i 



£2 = &,h + 



N h 



£ £ . 3^2,i£,2,h H3,i 

&-&,* + — flr— + ^r 



(29) 
(30) 



C4 =^ + t— + m + w (31) 



h ly h 

-2 



- _ - W^l2,iU,h . (10M3,i + 15jU2,i)6,h 
?5 — t,5,h H = r 



N h 



m 



1 (10/^2,1^3,1 + 5/l 4 ,»)$2,fe AfSji 



(32) 



etc. Clearly, the effort and complexity increase at each order. 
Identical results are obtained by the composition of gener- 
ating functions in equation (23). The general term can be 
obtained from the generating function K(t) for moments of 
total counts. Using the relation M(t) = G(t + 1), the com- 
position of probability generating functions in equation (23 
is also a composition of moment generating functions, 



K{t) = log[M (t)] = log[G(i + 1)] 

= \Qg{g h \9i{t + 1)]} = log{g h [Mi(t)]} 
= ]og{M h [Mi(t) - 1]} = K h [Mi{t) - 1], 



(33) 



from which it is clear that £*, continues to depend to all or- 
ders on the connected moments £ k ,h of halo number as the 
coefficients in Kh and the raw moments p,k,i of the halo occu- 
pation distribution as the coefficients in Mi . The generating 
function in equation (33) and the expressions for moments of 
counts in equations (29)-(32) plus extension to higher orders 
constitute the main result of the point cluster model. The 
point cluster results are independent of the internal details 
of halo profiles or concentrations. The general expression for 
£k contains contributions from occupation number moments 
of order 1 through k and halo correlations of order 1 through 
k; in £5, the first term arises from five objects in five separate 
haloes, the last from occupancy five in a single halo, while 
other terms represent four haloes with occupancies (2,1,1,1); 
three haloes with occupancies (3,1,1) and (2,2,1); and two 
haloes with occupancies (3,2) and (4,1). The numerical fac- 
tors represent the number of equivalent halo assignments. 
The sum of the combinatoric factors of a given £ Ui h in the 
expression for £*, are known as Stirling numbers of the sec- 
ond kind, S(n,k), the number of ways of putting n distin- 
guishable objects into k cells with no cells empty (Scherrer 
& Bertschinger 1991). Here they are produced from a gener- 
ating function in a manner such that any term desired can 
be easily produced by an algebraic manipulator. 

Some special cases are useful to consider. For single- 
element clusters, Ni — 1 with probability 1, the occupa- 
tion moments are pi = 1 and p, k — for k ^ 2, galax- 
ies are haloes and galaxy correlations are halo correlations, 
£fc = £k,h- For Poisson occupation number, the occupation 
moments are all p, kti = 1, and the halo model expressions re- 
produce the discreteness terms of equations (4)-(7). This is 
the locally Poisson realization of a distribution with spatially 
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varying n(r). For uncorrelated cluster positions, irreducible 
moments arise only from objects in the same halo; in this 
case the halo number Nh has a Poisson distribution, and the 
single-halo contribution to the count moment, 









(34) 



is often called the Poisson term. The composition of gener- 
ating functions for a Poisson halo distribution was studied 
by Sheth (1995a,b). The point cluster model varies from 
the halo clustering limit to the Poisson limit as a func- 
tion of scale. Typically, the two-point function behaves as 
£(r) ~ r~ 7 with 7 ~ 2, and higher order correlations scale 
hierarchically, as £& = Sk£ with nearly constant Sk- 
Since N grows as -R 3 , the dominant contribution to £,k on 
large scales then comes from the halo correlation, £*, « £k,h, 
but on scales where iV£ < 1, the point cluster model gives 
the one-halo term in equation (34). In this regime total 
number count moments have hierarchical correlations, with 

Sk = flk,i/fJ-2 

Many common statistical models are constructed start- 
ing with a Poisson halo number distribution, so that £k,h — 
for k ^ 2, and equation (34) holds exactly. If the occu- 
pancy distribution is also Poisson, Jlk,i = 1, then Sk — 1 
for all k (Si — S2 — 1 always), saturating constraints 
S2 m S2n ^ S^+n arising from the Schwarz inequality; this is 
a realization of the minimal hierarchical model of Fry (1985). 
Other examples of compound Poisson distributions include 
the negative binomial distribution, which is the composition 
of a Poisson cluster distribution with a logarithmic occupa- 
tion distribution (Sheth 1995b), and the thermodynamic or 
quasi-cquilibrium distribution of Saslaw & Hamilton (1984), 
which is the composition of a Poisson cluster distribution 
with a Borel occupation distribution (Saslaw 1989; Sheth &: 
Saslaw 1994; Sheth 1995a). 

There is one generalization that is also useful, where 
the total number of objects is the sum of contributions from 
two independent populations, N — N c + JVj, such as the sum 
of a strongly clustered population plus a weakly clustered 
"background" (cf. Soneira & Peebles 1977). In this case the 
cumulant moments £* are simply additive, 



N k ^k^N^ c , k + N^b,k. 



(35) 



If the background contributes to the total count, N — N c + 
Nb, but not to higher order moments, we have 



e* 



jv; 



(N b + N c ) k 



Sc.fc — Jc 4c 



(36) 



where f c is the fraction of clustered objects. Although the 
correlations are diluted, this says that the amplitudes for 
k > 3 are amplified, Sk = S c ,k/fc~ 2 - 



4 NUMERICAL RESULTS 

In this section we compare the model to statistics of galaxies 
and haloes identified within the setting of a single numerical 
simulation. The sample we use is the same one analyzed in 
Colombi, Chodorowski & Teyssier (2007, hereafter CCT), 
where many more details can be found. The simulation is 
performed with the adaptive mesh refinement (AMR) code 



RAMSES (Teyssier 2002), assuming a standard ACDM cos- 
mology with n m = 0.3, Q.a = 0.7, Ho = lOO/tkms -1 Mpc" 1 
with h = 0.7, and a normalization as — 0.93, where as is the 
root mean square initial density fluctuations in a sphere of 
radius 8ft -1 Mpc extrapolated linearly to the present time. 
The simulation contains 512 3 dark matter particles on the 
AMR grid, initially regular of size 512 3 , in a periodic cube 
of size Lbox = 200 h^ 1 Mpc; the mass of a single particle is 
then 7.09 x 10 9 Mq . Additional refinement is allowed during 
runtime: cells containing more than Namr = 40 particles are 
divided using the standard AMR technique with a maximum 
of 7 levels of refinement. 

A halo catalog, Eh, and a "galaxy" (subhalo) catalog, 
Eh, are extracted from the final state of the simulation us- 
ing the publically available software adaptaHDP (Auber, Pi- 
chon & Colombi 2004); details of the procedure can again be 
found in CCT. We use the number of dark matter substruc- 
tures in each halo detected by adaptaHOP as a proxy for the 
galaxy distribution. A halo can contain one or more galax- 
ies: a single component halo hosts one galaxy (or is its own 
substructure), and an TV-component halo hosts N galaxies. 
The substructure distribution differs somewhat from that of 
galaxies (see the discussions in CCT and Weinberg, Colombi, 
Dave & Katz 2008), but it provides a discrete number count 
distribution that is useful to test how the behavior of the 
discrete halo model differs from that of the continuous mass 
distribution. 

Figure 1 shows the distributions of halo mass f(m) 
and of occupation number Pn for the full halo catalog. 
The range of masses covers almost four decades; the largest 
halo contains 53 substructures. Moments of these distri- 
butions give the occupation moments JXk that appear in 
equations (29)-(32) and the point halo hierarchical am- 
plitudes Sk — flk/fi-2^ 1 - F° r comparison, smooth lines in 
the figure show the Press-Schechter (solid line) and Sheth- 
Tormen (dotted line) mass functions, plotted for S c = 1.50. 
The Press-Schechter and Sheth-Tormen mass distributions 
provide a good representation of the mass function for 
M > 10 12 Mq, rising with mass a little more weakly than 
1/M towards small masses and with an exponential cutoff 
at large mass. The number distribution behaves as a power 
law, Pjv ~ 1/N P , with p in the range 2.0-2.4. The subclump 
finder adaptaHOP identifies haloes as connected regions with 
density contrast larger than S > 80 employing a standard 
SPH softening of the particle distribution with Asph = 64 
neighbors (see, e.g., Monaghan 1992). This, along with the 
mass resolution of the structures resolved by RAMSES, con- 
trolled by the value of JVamr, leads to the rather soft small- 
M cutoff on the halo mass function in Figure 1 and also the 
low value of <5 C . 

The full samples Eh and E g contain 50234 haloes and 
64316 "galaxies", respectively. Most haloes have a single 
component; the average number of substructures per halo 
is Ni = 1.28. From these two parent catalogs, we apply var- 
ious mass thresholdings to extract subsamples from Eh and 
Eg that we denote _Eh(M m i n , M max ) where M m i n and M max , 
given in solar masses, correspond to minimum and maxi- 
mum mass thresholds of the host haloes respectively. We 
use these subcatalogs to test the variation of halo clustering 
with mass. The different realizations break down as follows: 

(i) The full sample separated into "light" and "mas- 
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Figure 1. Distribution of halo masses f(M) (histogram; bottom 
scale) as a function of M and substructure occupation number 
probability Ppf (shaded histogram; top scale) as a function of 
N. The solid and dashed curves show the Press-Schechter and 
Sheth-Tormen mass functions. 



Table 1. Details of the catalogs extracted from RAMSES. The two 
first columns give the halo mass range, Af m ; n < M h < M max (in 
units of Mq); the third and fourth columns give the number of 
objects in the halo and substructure catalogs; and the fifth and 
sixth columns give the average number of substructures and dark 
matter particles per halo. 
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sive" halo subsamples, Eh(0, oo) = Eh = Eh(0, 10 14 ) + 
£ti(10 14 , oo), and the substructure counterparts. 

(ii) A catalog of haloes with masses larger than 5 x 
10 11 Mq, which avoids the strongest rolloff at small mass, 
separated likewise into two subsamples, -Eh (5 x 10 11 ,oo) = 
E h (5 x 10 11 , 10 14 ) + E h (10 14 , oo), and the substructure coun- 
terparts; 

(iii) A catalog of haloes with masses larger than 4 x 
W 12 Mq, which avoids essentially all of the rolloff at small 
mass, separated likewise into two subsamples, -Eh (4 x 
10 12 ,oo) ee E h (4 x 10 12 ,10 14 ) + -E h (10 14 ,oo), and the sub- 
structure counterparts. 

Table 1 summarizes subcatalog information. 

From the distribution of mass and the catalogs of haloes 
and subhaloes in the simulation we compute correlation 
statistics |j for k = 2-5. Figure 2 shows the variance, or 
volume-averaged two-point correlation function £2, evalu- 
ated for spherical volumes of radius R as a function of 
R, for dark matter, or mass (solid line), and for haloes 
(long-dashed line). Dotted lines show the predictions of the 
point cluster model: the upper line for only the mass in 



haloes, and the lower line including mass not contained in 
haloes as an unclustered background, as in eq. (36), with 
f c = y] Mh /Mtot = 0.36. Finally, the short-dashed line in- 
cludes two additional aspects from the full halo model, a 
modest relative bias factor b — 1.22 between mass and haloes 
on large scales, and the effects of resolved haloes, detailed 
in Section 5 below. 

Panels in Figure 3 show the second moment evalu- 
ated for substructure "galaxies" £2,9 (solid lines) and for 
haloes, &,h (long-dashed lines) for the four inclusive cata- 
logs: haloes of all masses, haloes with M > 5 x 10 11 Mq, 
with M > 4 x 10 12 Mo, and with M > 10 14 Mr*, as iden- 



*©> 



©> 



tilled in the caption. Dotted curves show the predictions of 
the simple point cluster model. Although it has shortcom- 
ings in detail, on scales of a few Mpc the point cluster model 
with no adjustments reproduces the trends with scale and 
from catalog to catalog to within a factor of two or so over 
four decades of correlation strength. Dashed curves show a 
quantitative improvement with a very modest adjustment of 
parameters, relative bias factors of 1.2 or 1.3 and occupation 
moments adjusted by a factor of 2, as given in the middle 
columns in Table 2. For radius smaller than 1 h~ 1 Mpc finite 
halo size starts to become important. The mass threshold 
M > 5 x 10 11 Mq removes only haloes containing a sin- 
gle substructure (the smallest halo containing two substruc- 
tures has a mass 8 x 10 11 Mq), and so affects only the mean 
N = (N\ but none of the higher factorial moments / 7V' fe ' V 
just as for an unclustered background population as in equa- 
tion (36). Thus, in the regime where the normalized moment 
is large, |jb S> 1, it is simply rescaled, |jj./|jb = (N/N') k . 
This is apparent for the data plotted in Figure 2, where 
for the smallest cells £2 for the E g (5 x 10 11 ) subsample is 
larger than that for the full E s sample by a factor 1.249, 
very close to the number ratio (64316/57564) 2 = 1.246. The 
next mass threshold, M > 4 x 10 12 Mq, removes doubly 
and also triply occupied haloes (the smallest halo contain- 
ing three substructures has a mass 1.8 x 10 12 Mq), and so 
this threshold changes the shape of £2 and £3. 

Panels in Figure 4 show the hierarchical amplitude S3 
for the four subhalo catalogs (solid lines), and for the cor- 
responding halo catalogs (long-dashed lines). Finite volume 
limitations are apparent at large scales, and the Eh(10 14 ) 
sample is not large enough for a reliable third moment on 
almost any scale. Dotted lines show the naive point clus- 
ter model. Again, the first mass cut M > 5 x 10 Mq 
removes only haloes containing a single substructure from 
the full catalog, changing only the mean count; and on 
the smallest scales the expected scaling S3/S3 = N'/N = 
57563/64316 = 0.895 is again satisfied. 

Figure 5 shows the amplitudes S3, Si, and S5 for for 
dark matter (solid lines) and for haloes (long-dashed lines) . 
Figure 6 shows the St for substructures (solid lines) and for 
haloes (long-dashed lines), for the entire Eh halo sample. 
The naive point cluster model agrees with the simulations 
qualitatively but not quantitatively. One possible explana- 
tion is that halo occupation is correlated with environment, 
and a modest adjustment of the point cluster parameters 
gives a good fit. Table 2 shows the naive point cluster model 
result using occupation probabilities pn and the halo mass 
function n(M) from the simulations, and also the result of 
adjusting fit parameters. In the point cluster model, the pa- 
rameters are factorial moments /Zfc = ( N<- k > \ / N k for galax- 
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Figure 2. Variance £2 (?") for mass (solid line) and for halo number 
(long-dashed line) measured for spheres of radius r in numerical 
simulations compared with the halo model. Dotted lines show the 
point cluster model of equation (29); the upper line shows only 
the contribution of particles in haloes, while the lower dotted line 
includes particles not in haloes as an unclustered background. 
The short-dashed line includes a bias b = 1.25 between mass and 
haloes on large scales and the effects of resolved haloes on small 
scales. 
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Figure 3. Cell count variance £2(7*) in spheres of radius r ob- 
tained from numerical simulations compared with the halo model. 
Solid lines show £2 measured for galaxies (substructures) and 
long-dashed lines show £ 2j j, measured for haloes identified in the 
simulations for the four subcatalogs. Results for all haloes are 
presented in panel (a), M > 5 X 10 11 in panel (b), 4 X 10 12 in 
panel (c), and 10 14 Mq in panel (d). The dotted lines show the 
point cluster model of equation (29), and the short-dashed lines 
include bias and resolved haloes. 



ies and fik = ( M k \ /M k for mass, computed for the haloes 
identified in the simulation. The quantity identified as "6" 
is the large-scale relative bias between galaxies and haloes, 



Figure 4. Hierarchical amplitudes Ss(r) from numerical simula- 
tions for substructure "galaxy" catalogs (solid lines) and for the 
corresponding haloes (long-dashed lines) The dotted lines show 
the naive point cluster model, and the short-dashed lines include 
bias and resolved haloes. Panels show the four different subcata- 
logs, as in Figure 3. 
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Figure 5. Hierarchical amplitudes Sh(r) for mass density from 
numerical simulations for k = 3, 4, and 5 (bottom to top). Solid 
lines show Sk measured for mass; long-dashed lines show S^ for 
haloes. Dotted lines show the naive point cluster model expres- 
sions in equations (30-32) and same results but adjusted for a 
weakly clustered background (not contained in haloes). The short- 
dashed lines also include a bias factor on large scales and the ef- 
fects of resolved haloes on small scales, as detailed in Section 5 
below. 



5 THE FULL HALO MODEL 

To extend our understanding we turn to the context of 
the recently developed phenomenological halo model (Seljak 
2000; Ma & Fry 2000b; Peacock & Smith 2000; Scoccimarro, 
Sheth, Hui, & Jain 2001; Cooray & Sheth 2002), in which 
clustering on small scales is derived from the mass function, 
profiles, and clustering properties of dark matter haloes. Nu- 
merical simulations have suggested haloes have a universal 
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Table 2. Correlation parameters for "galaxies" and mass. The first five columns show values in the simple point cluster model, in which 
all haloes are identical, using the simulation halo distributions pjv and f(M); the next five columns show fit values (f); and the last five 
columns show values computed in the detailed halo model (m). 
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°5 


6/6- 


A? 


cm 
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cm 


S™ 


dark matter (mass) 


1 


4.00 


7.18 


93.0 


1460 


1.25 


1.96 


7.2 


93 
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1.71 


34.3 
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2450 
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1 


1.35 
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3360 


1.22 
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m > 5 X 10 11 Mq 
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300 
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6.95 


5.24 
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m > 4 X 10 12 Mq 


1 


2.10 


4.00 


26.8 


212 


1.34 


1.27 


2.3 


10 


68 


1.30 


3.03 


3.27 


22.4 
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m>lx 10 14 Mq 
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1.50 


2.93 


6.59 
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where W(x) — 3(sina; — :rcos:r)/:r 3 is the Fourier window 
corresponding to a real-space top-hat window function, and 
m = 47tp_R 3 /3. In terms of v, the density of haloes of mass 
M is then written 



"/(") 



(39) 



dn d(ln v) p 
Am d(lnm) m 2 

The Gaussian distribution function f(y) of Press & 
Schechter (1974, PS) and the refinement of Sheth & Tor- 
men (1999, ST) have 

!//(>/) =2A[l + (qv 2 T P ] (£) 1/2 e"^ 2/2 . (40) 



R (h-'Mpc) 



The normalization A is chosen so that J Av f(y) — 1 and is 
independent of q. The Press-Schechter function has q — 1, 
p — 0, A — |; the Sheth-Tormen form has q = 0.707, p — 
0.3, A « 0.32218. 



Figure 6. Hierarchical amplitudes Sk{r) for galaxy number count 
measured in from numerical simulations for k = 3, 4, and 5 (bot- 
tom to top). Solid lines show galaxy (substructure) Sfc, and long- 
dashed lines show halo Sk,h- Dotted lines show the point cluster 
model expressions in equations (30-32). 



density profile (Navarro, Frenk, & White 1997; Moore et al. 
1999; Navarro et al. 2004), 



P(r) 

P 



My) 

I s 



(37) 



where the scale r 3 and amplitude A are functions of the halo 
mass. In particular, r s is related to the virial radius 7-200, 
within which the average density is 200 times the mean, by 
a "concentration parameter" c(m), r s = r2oo/c; this then 
also determines the amplitude A. For a large cluster, say 
m = 10 15 Mq, the virial radius is about 3 /i _1 Mpc; and with 
c ~ 6 the scale radius is r s ~ 500/i~ 1 kpc. Thus, at roughly 
Mpc and smaller scales we begin to resolve clusters, and we 
expect to have to replace the point cluster model with the 
full halo model. 

The halo mass function is conveniently written as a 
function of the dimensionless overdensity u = S c /a(m), 
where S c is the threshold overdensity that leads to a col- 
lapsed halo, often <5 C = 1.68, and cr 2 (R) is the mean square 
mass fluctuation within a sphere of radius R evaluated for 
the linearly evolved input power spectrum. Specifically, 



a (m) 



d J fc 
(2tt) 3 



P(k)W 2 (kR), 



(38) 



5.1 One-Halo Term in the Point Cluster Limit 

In the point cluster limit of the full halo model, the sum 
over haloes in equation (16) and the resulting composition 
of generating functions in equation (23) remain true, but 
the calculation now includes an average over the distribu- 
tion of halo masses as well as over halo occupation and halo 
count, both of which now differ with halo mass. The result- 
ing order-fc connected correlation function is again a sum of 
contributions from a single halo to k different haloes, just as 
in equations (29)-(32), with the same coefficients. For the 
one-halo term of the order-fc moment in the full halo model, 
the average over all haloes includes an average over the halo 
mass function dn/dm, 

( N l k] ) I dm (dn/dm) V ( N [k] (m) ) 






(Ni) [J Am (dn/dm) V ( N{m) )] 



(41) 



where N(m) — Ni(m) is the occupancy of a halo of mass m. 
The factor in brackets in the denominator is 



An 



Am — V { N(m) )=n h V(N l )= N h N t 



N, 



(42) 



where fih is the number density of all haloes, Nh — nnV 
is the mean number of haloes in V, and Ni is the average 
occupation over haloes of all masses; and the numerator is 



dm — V(N lk] (m)) 
Am x ' 



nuV(N^) 



N h N k fi k , 



(43) 



The one-halo term of the full halo model thus produces 
the same result as the previous point cluster result, ^ h = 
Mfe.i/iV^ -1 , with occupation moment 
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(j-h,i = 



J dm (dn/dm) ( N [k] (m) ) / J dm (dn/dm) 
[J dm (dn/dm) ( N(m) ) / Jdm (dn/dm)] * 



(44) 



We can compute occupation moments flk for mass from 
first principles by taking number N to be proportional to 
mass (the number of hydrogen atoms or dark matter parti- 
cles), N on m and N*-' oc m fc (with no discreteness terms). 
From (44), these are 



/<fc 



[ J dm (dn/dm) m k ] [J dm (dn/dm) 



[J dm (dn/dm) m ] k 
with corresponding hierarchical amplitudes 



Sk = 



Mfc 

at 1 



[Jduf(u)m k - 1 ][fduf( V )]' 

[Jduf^m]"- 1 



(45) 



(46) 



over the range of scales where the one-halo term dominates 
but haloes are not resolved. In this case, results are deter- 
mined entirely by the mass function, which in turn is related 
to the primordial power spectrum. For a power-law power 
spectrum, with v — (m/mi)"" 1 "™" 6 , and with the Press- 
Schechter and Sheth-Tormen forms of the mass functions, 
the integrals can be done analytically, giving 



S k = 



I(k) [1(1)} 
[1(2)]*-: 



where 



i(k) = r 



3(fc 



3 + n 



+ 2- 



3(fc-l) 



3 + n 



(47) 



(48) 



independent of q. For a Poisson cluster distribution (on 
small scales cluster correlations are unimportant) and for 
the Press-Schechter mass function, this expression was also 
obtained by Sheth (1996b). For the Press-Schechter mass 
function, which has p — 0, and for spectral index n = this 
gives the particularly simple result Sk — (2fc — 3)!!. Results 
for power-law spectra are shown in Figure 7, together with 
results from numerical simulations by Colombi, Bouchet, & 
Hernquist (1996) (plotted are the values of Sk measured 
at £ = 100, but values at £ = 10 or £ = 1 differ by less 
than the error bars). The Sheth-Tormen mass function ap- 
pears to agree poorly with the numerical results; this is one 
instance where the observed behavior seems to prefer the 
Press-Schechter form, at least for n not too negative. How- 
ever, the Sheth-Tormen function is relatively more weighted 
towards smaller masses, and in numerical simulations there 
is always a smallest mass that can be considered. Thus, we 
examine the results of a small mass cutoff in the integral, of 
10~ 4 and 10~ 2 in units of the mass m\ at which v(m\) = 1. 
A 10~ 4 cutoff mass has little effect on PS but is significant 
for ST, and a 0.01 cutoff has a significant effect on both. In 
the simulations, the ratio of the particle mass to the non- 
linear mass is typically in the range 0.001-0.01, and the ST 
mass function with a moderate low-mass cutoff does agree 
with the simulations results, at least for —1 <n< 1. Asn 
becomes more negative, all the halo model curves rise much 
more rapidly than the trend seen in the simulation results. 
This may reflect an increasing difficulty in simulating nega- 
tive values of n (cf. Jain & Bertschinger 1998) . 




Figure 7. Amplitudes Sk for k = 3, 4, 5 (bottom to top), as 
a function of spectral index n for power-law spectra. For each k 
curves are two sets of curves; the lower set (blue in color) shows 
results for the Press-Schechter mass function, and the upper set 
(red in color) for Sheth— Tormen, computed from eq. (46). Solid 
lines show Sk integrating over all masses; long-dashed curves have 
a lower mass cutoff m > 10 — 4 ; and short-dashed curves have 
m > 10~ 2 , in units of the mass ra\ at which v(m\) = 1. The 
dot-dashed (green) curves show the predictions of hyperextended 
perturbation theory (Scoccimarro & Frieman 1999). Symbols with 
error bars show results from numerical simulations (Colombi, 
Bouchet, & Hernquist 1996). 



For statistics of galaxy number counts, we must average 
over moments of halo occupation number, 

s _ [J du m ( N [k] (m) ) /m] [/ dv f(v) (N(m)) /m] ^ 
[fdvf(v)(NW(m))/m] k - 1 

(49) 

(the factor 1/m remains from the PS or ST halo mass func- 
tion). In simulations, in general it is found the mean num- 
ber of galaxies / N(m) \ grows more slowly than linearly 
in mass. Models have included a power-law, / N(m) \ = 
(m/mi) , with /3 < 1 and perhaps with a minimum mass 
cutoff mo; a broken two-power-law model (Berlind & Wein- 
berg 2002); and a similarly behaved but smoothly interpo- 
lated function (Berlind et al. 2003). Substructures or subhalo 
occupation numbers exhibit a similar behavior, but perhaps 
with /3 — > 1 at high mass (Kravtsov et al. 2004) . Higher order 
correlations also require higher order moments of the halo 
occupation distribution, which are typically sub-Poisson at 
small N, with / N(N - l)(m) ) < ( N(m) ) . From semian- 
alytic galaxy formation considerations, Scoccimarro, Sheth, 
Hui, & Jain (2001) extend to higher orders by assuming a 
binomial distribution, also used by Kravtsov et al. (2004). 
However, a representation of the galaxy number count dis- 
tribution as a central galaxy plus a Poisson distribution of 
satellites (Berlind et al. 2003; Kravtsov et al. 2004; Zheng et. 
al 2005; Zheng & Weinberg 2007; Zheng et al. 2007, 2009) is 
becoming increasingly popular. In this representation, both 
central and satellite distributions are characterized entirely 
by their means \N C ) — N c and \N S ) — N s . Since iV c takes 
on only the values and 1, so that for any positive power 



Cell Count Moments in the Halo Model 9 



/ N£ \ — N c , and since for a Poisson satellite distribution 
\Ng ) — N s , the factorial moments of halo occupancy are 

( N lk] ) = iV* +k N c N*-\ (50) 

The central object can be modeled as a sharp or smoothed 
step function (Berlind et al. 2003), and Zheng et. al (2005) 
present expressions for moments with parameters extracted 
from simulations. A main import of all models is that mo- 
ments of occupation number grow more slowly than linearly 
with mass, a behavior that we model as the simpler form 
N(rn) ~ mr with j3 < 1. With the Sheth-Tormen mass func- 
tion and with no small-mass cutoff, Sk for number counts is 
again a ratio of T-functions, 



S, 



im \h{i)\ k - 



where now 

3/3(fc-l) 



iff 



(3 + n) 



1 

+ 2. 



+ 2 p r 



3/3(fc - 1) 



(3 + ' 



2- p . 



(51) 



(52) 



We can use the Poisson model to obtain the occupancy 
probability distribution averaged over all haloes. For a Pois- 
son distribution with mean /x(m), the probability pn for a 
halo of mass m to contain N galaxies (or N satellite galaxies) 
is pn — fi N e -M /JV!. Averaged over the power-law portion 
of the Press-Schechter mass function dn/dm ~ u/rn 2 , with 
the integral cut off by the Poisson exponential e~ M before 
the exponential cutoff in the mass function is reached, the 
probability pn of N objects in any halo scales as 

[N+(3 + n)/6/3-l-l//3]\ 



PN CC 



AT! 



where v ~ m (3+n)//6 and \x,(m) 
this behaves as a power law, 



J 



PN 



N~ 



W 



(3 + n) 
6/3 ' 



(53) 
As N becomes large, 

(54) 



For n ~ —2 and /3 < 1 the exponent is near r — —2, a good 
approximation to the distribution plotted in Figure 1. For 
Sheth-Tormen the power is shifted by 2p(3 + n)/6/3, or by 
about 0.1. 



5.2 Resolved Haloes 

For small volumes we can no longer take haloes as point ob- 
jects, but must take into account the distribution of objects 
within a halo. In the full halo model, the one-halo contribu- 
tion to the fc-point function ££ for mass is a convolution of 
halo profiles (Ma & Fry 2000b), 






/ dm (dn/dm) m k J A 3 r'u{y' 1 ) ■ ■ ■ u{y' k ) 
\J dm (dn/dm) m J d 3 r' u(y')] 



(55) 



where the position r' of the halo centre runs over all space, 
y[ = \n — r'\/r 3 , and the scaled halo profile u{r) is nor- 
malized to unit integral. From equation (55), the volume- 
averaged correlation is then 






J dm (dn/dm) m fe / dV [F{r')] k 
[J dm (dn/dm) m J d 3 r' F(r')] k '' 



(56) 



where F(r') is the portion of the total volume of a halo 
centred at r' that lies within V, 



F(r') = 



(57) 



Note that the integrand is a function of r/r 3 , and since the 
scale radius r a depends on mass, the from factor F is in 
general also a function of halo mass. From equations (56) 
and (57) we can recover the point cluster model: if a volume 
is much larger than a halo size, R S> r s for all haloes, then 
F(r') is very small unless the halo itself is within V, in which 
case the integral then contains the entire halo contents. In 
this limit and with unit normalization, F — > 1 for r' in V 
and F — > for r' outside V . Then, the integral over r' is 
just a factor of V, and we recover the point cluster model. 

For resolved haloes in moments of discrete galaxy counts 
we consider first the second count moment £2- Let a halo 
contain iV objects, and let N' be the number of these objects 
that are contained within V . Then N' = ^2 Ni, where either 
Ni — 1 with probability pi if object i is counted or Ni — 
if object i is not, and the second moment is 



JV JV 



(N'(N'-1)) = (J2N 1 (Y / N 3 -1) 
x »=i j=i 



N N 

1-2 






(58) 



i=l i=l 

But since Ni takes on the values or 1, Nf — Ni and the last 
two terms cancel, leaving the sum only over distinct objects 



( N'{N' - 1) > = ( £ NiNj ) = (N(N-l))p 2 



(59) 



If object positions within a halo are uncorrelated, the prob- 
ability p that an object within a given halo is located within 
the volume V is just the fraction F of the halo that is within 
V, form factor in equation (57), the same for all objects and 
independent of the halo occupation JV, 

(N'(N'-l)) = (7V(A-1))(F 2 ). (60) 

This agrees with the usual practice, to distribute the average 
pair count ( N(N — 1) V weighted by the square of the halo 
profile form factor (F 2 \, 



/'2 



/ dm (dn/dm) ( N(N - 1) ) (F 2 ) / n h 



(61) 



(62) 



[ / dm (dn/dm) ( N(m) ) / n h ] 
where the volume-averaged form factor is 

( Fk ) = vJ ^[W] ■ 

In the position space formulation symmetry over all par- 
ticles is manifestly maintained in the form-factor integrals, 
without need to introduce the approximation W12 ~ W1W2. 
The form factor F does not appear in iV in the denominator 
of equation (61), since, as can be easily seen by changing the 
order of integration, \F) = 1. The calculation for a halo 
occupation distribution consisting of a central object plus 
N s — N — 1 satellites yields 

( N' [2] ) = (n1 2] )(f 2 )+2(N s )(fF c ), (63) 

where F c — 1 for r < R and vanishes otherwise. Extending 
to general fc, we obtain 

( N' [k] ) = ( N [k] ) ( F k ) (64) 
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with no central object, or 

< N' [k] ) = < JVJ* 1 ) ( F k ) +fc ( TV?" 11 ) ( F k -' F c ), (65) 

with a central object. The last term could contain a factor 
\N C \ if this is not 1. The form- factor-corrected halo occu- 
pation moment is then 



Mfc 



/ dm {An /Am) ( N [k] (m) ) ( F k ) / n h 
[ J Am {An/ Am) ( N{m) ) / nh ] , 



(66) 



modified as in equation (65) for a central object. 

For moments of dark matter mass, a reasonably good 
representation of the numerical results is obtained using 
the NFW profile, but for substructures this is not the case. 
The substructure profile was seen in Diemand et al. (2004) 
to follow roughly an isothermal profile, and we have stud- 
ied using the isothermal sphere profile also. The measure- 
ments of Diemand et al. (2004) and our own do not pro- 
vide enough statistics to infer a mass dependence of the 
concentration parameter, and so we use a constant value 
c = 10 that gives reasonable results on small scales. Fig- 
ure 8 shows the volume-averaged form factor for fc = 2-5 
for NFW haloes (solid lines) and for the isothermal sphere 
profile (long-dashed lines), both with c = 10. Curves are 
plotted as a function of Y — R/r s , where r s = r2oo/c. As 
expected, the form factor goes to 1 at large scale and falls 
rapidly for small R, where only a small fraction of a halo is 
sampled. Note that ( F k ) < ( F" ) if k < n. The integral 
converges to 1 on large scales, the point cluster regime, but 
falls rapidly for Y < c. In equation (66), for fixed R, this 
factor decreases rapidly for increasing mass. 

Figure 9 shows the form-factor corrected, one-halo Sk — 
P-k/p-^" 1 , normalized by its value in the point-cluster limit, 
as a function of R, for k — 3, 4, and 5 (bottom to top; dif- 
ferent orders k offset for clarity). On small scales, smaller 
than a few Mpc, this shows the effect of resolved haloes. 
The result depends on both halo profile and on the distri- 
bution function: solid lines show NFW haloes averaged over 
mass; long-dashed lines show the same haloes averaged over 
number; short-dashed lines show the isothermal profile with 
c = 10; and dotted lines show isothermal profiles with con- 
centration c{m) as for NFW. On large scales halo profile 
shape has no effect, but on small scales the differences for 
different profiles and weightings are substantial. 

From the halo model can extract small-scale behaviors 
of the correlations £&, which can be different for galaxies 
and for mass. The concentration parameter plays a critical 
role in the result. For scale invariant spectra (c.f. Davis & 
Peebles 1977) we expect c{m) ~ M~ a , with a = (3+n)/6 (in 
Ma & Fry 2000a, this parameter is 0). For ACDM, over our 
relatively small range in mass we also take the concentration 
parameter to scale as a power of mass, a ~ 0.11 or 0.12 
(Bullock et al. 2001; Zhao et al. 2003), corresponding to an 
effective n ~ —2.3. As above, let the number of objects in 
a halo grow with mass as / N^ ' \ ~ m k ^ (for statistics of 

mass f3 = 1). Finally, let An/ Am ~ v v /m 2 as m — >• 0, where 
v ~ m (3+n)/6 and p' = 1 for PS and p' = 1 - 2p = 0.4 
for ST (in Ma & Fry 2000a, this parameter is a). Then, 
ignoring the exponential factor in An/ Am on small scales 
where v is small and changing integration variable from m 
to Y — R/r s = cR/r 200 in equation (66), we see that p,k 
scales as 




Figure 8. Volume-averaged form factor ( F k \ for k = 2-5 (upper 
left to lower right), as a function of Y = R/r s . Solid lines shows 
result for NFW profile; long-dashed lines show isothermal profile, 
both with c = 10. 



^k 



[* 



l/( Q +l/3)]['=/ 3 +P'(3+™)/6-l] 



(67) 



and the fc-point function £& = /ik/N k x scales as R 7fc , with 



7fc 



l + 3a 



[3(fe-l)a + fe(l-/9)] 



(3 + n)p' 
2(1 + 3a) 



(68) 



= 3(5-2/3 + n) 6(l^_(3 + ny 

v ' 5+n 5+n 5+n v/ 

independent of the shape of the halo profile. For /3 = 1 this 
is the same as the result obtained by Ma & Fry (2000a) 
(beware a change of notation) for mass, and for fi — 1 — e 
is the result obtained by Scoccimarro, Sheth, Hui, & Jain 
(2001) for galaxy number. This is of the hierarchical form 
only for p = 0, which is not true for either of the PS or ST 
mass functions, and for /3 = 1. Departures from hierarchical 
scaling in the small- R behavior of Sk grow with k, 



where 



A7 



(3 + n)p' - 6(1 - 13) _ (3 + n)p' - 6(1 - g) 



2(1 + 3a) 



5 + n 



(70) 



(71) 



(A7 « -0.26 for the choices p' = 0.4, /3 « 0.8, n « -2). The 
presence of ever higher powers of fa in Sk = ffc/?2 _1 em " 
phasizes any scaling defects in £2- An interesting alternative 
normalization is 



S' k = 



& 



- (k -l)/(k-2) 
Sfe-1 



(fe_l)/(fc_2) 
°k-l 



(72) 



Departures from scaling in S' k decrease with k for fc ^ 3, as 
S' k ~ i? A ^( fc " 2 ) (73) 

for the same A7 given in equation (71). 

5.3 Multiple-Halo Terms 

Terms that involve objects in multiple haloes also depend on 
correlations among haloes. In the perturbative regime, halo 
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Figure 9. Form-factor corrected one-halo Sk as a function of 
_R, normalized to its point cluster value, for fc = 3, 4, and 5 
(bottom to top; different k offset for clarity). Solid lines show S^ 
for NFW profiles weighted by mass m fe ; long-dashed lines show 
NFW profile weighted by number N^'; short-dashed lines show 
isothermal profile with c = 10 weighted by number; and dotted 
lines show isothermal profile with c = c(m). 



correlations have bias factors that are functions of the halo 
masses, and higher order correlation functions also involve 
higher order bias parameters (Fry & Gaztahaga 1993); for 
instance, the halo three-point function is 



&,h = b(m 1 )b(m2)b(m 3 )^3,p(r 1 2,r 1 3,r23) 

+ b 2 fca,„(ria)$i,p(ris) + eye. (123)], 



(74) 



where the £ p are correlation statistics for the underlying (pri- 
mordial) density distribution. As a function of mass, the 
linear bias factor found for PS by Mo, Jing & White (1996, 
1997) and adapted for the ST halo mass function (Sheth & 
Tormen 1999; Casas-Miranda et al. 2003) is 



1 + 



qv 



1 2p 

- + 



S c [1 + (qv 2 )?] 



(75) 



with further refinements for small mass suggested by Jing 
(1999). Higher order functions also require higher order bias 
parameters (Mo, Jing & White 1997; Scoccimarro, Sheth, 
Hui, & Jain 2001), 



b 2 = 



21<5 C 



+ - 



qv 



2 4 
q v 



1 + 



2p 



, 2 W 



3qu + 



[1 + (r 

2p(qu 2 + 2p 



1) 



[1 + (qv 2 



(76) 



etc. Higher order bias terms 63, etc., vanish when integrated 
over the full halo mass function. Even with a low mass cutoff 
or with different mass or number weightings we expect that 
they remain generally small; and so we will drop them from 
now on (but see Angulo et al. 2008). 

We exhibit in detail the fc-halo and two-halo contribu- 
tions to £k in the full halo model. The fc-halo contribution 
to £k is 



N k e k h =n/ 



dm 






*>> 



(77) 



x / d 3 r[u(y[)---d ,i r' k u(y' k )£ kth (r' 1 , 



where N is as given in equation (42). Ignoring non-linear 
bias terms, so that in terms of the underlying density 
correlation £ k>p the halo correlation function is £k,h = 
b(mi) ■ ■ ■ b(rrik) £k,p, the volume-averaged correlation be- 



N' 



'tk h = i[[ dm >-^- & K) < N( mi ) ) 



x d 3 r' 1 F(r' 1 )---d 3 r' h F(r' k )Z klP (r'i, 



(78) 



,r h 



In the point cluster limit on large scales, for which F — 1 
for r' in V and F — for r outside V, this gives 



skh 
£,k,h 



5*6 



k,p 



uk $k,h, 



with an occupation-number weighted bias factor, 
/ dm (dn/dm) ( N(m) ) b(m) 



b = 



J dm (dn/dm) ( N(m) ) 



(79) 



(80) 



The halo correlation function £k,h is £,k,h = b h £k, P , with a 
bias factor weighted only by the halo mass distribution, 

J dm (dn/dm) b(m) 



b h = 



J dm (dn/dm) 



(81) 



Factors of mass or number weight greater contributions at 
higher masses, where b(m) takes on larger values, so in gen- 
eral b >bh\ galaxies are more strongly correlated than haloes 
on large scales, though only by a small amount. Ratios of 
integrals b/bh over the Sheth-Tormen mass function with 
/ A ^ oc m for mass and ( N \ oc nnr with /3 = 0.8 for num- 
ber are listed in Table 2. 

Similarly, we can write intermediate terms. The two- 
halo contribution to £& is a sum of terms of the form 

x / d 3 r[ dM [F(ri)] fcl [F(r' 2 )] k * b(m,)b(m 2 ) &(r 12 ), 



(82) 

where k = ki+k 2 . (If /ci = k 2 there is an additional symme- 
try factor of s — i because the partition and its complement 
are identical; the generating function gives all combinatoric 
factors automatically.) On large scales, where the halo size 
is insignificant and the form factors take the value F = 1 
over essentially the entire volume V, the full two-halo term 
is thus the sum over partitions 



J2h 



sk\ bk^bk2_ Mfci "fc 2 $2,fe 
kilk 2 l b h b h N k ~ 2 ' 



where bk is weighted by / N^ k ' V 

/ dm (dn/dm) ( N [k] (m) ) ( F k ) b(r 



bk = 



J dm (dn/dm) ( AM ) ( F k ) 



(83) 



(84) 



For moments of mass, the factors (N'- k '(m)) become m k . 
Weighted by different factors of number or mass, the bias 
parameters bk will in general be different from 6 = 61 defined 
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in equation (80); the lower mass limit for the integral also 
increases for higher order moments. On small scales, where 
haloes are resolved, the halo size, and thus the factor F(r'), 
depend also on halo mass: the mass and position integrals 
cannot be factored or simplified. However, since F ^ 1, the 
expression in equation (83) is an upper limit to the two-halo 
contribution, and even without the form factors the two- 
halo contribution is dominated by the one-halo term given 
in equation (56) as R — s- 0. Extension to other intermediate 
orders follows similar lines. 



6 DISCUSSION 

We have studied the behavior of cell count moments, in- 
cluding the variance £2 and the hierarchal amplitudes Sk for 
k = 3, 4, and 5, in the context of the halo model, and we 
have compared the model with results of numerical simula- 
tions for statistics mass and of galaxy (substructure) num- 
ber counts identified in the same simulation. Expressions 
(29)-(32) constitute the halo model predictions for the two-, 
three-, four-, and five-point functions; a composition of gen- 
crating functions for the halo number and halo occupancy 
distributions, as presented in equation (33), produces auto- 
matically the halo model result at general order, including 
all terms and combinatoric factors. The naive, point-cluster 
form of the model with identical haloes is easily generalised 
to include averages over a distribution of halo masses and 
over positions within resolved haloes. The general form of 
the naive point cluster model results continues to hold, with 
the addition of a modest bias, of a factor of two or less, 
on large scales, and form factors that reflect shapes of re- 
solved haloes on small scales. With these components, the 
halo model is able to reproduce in quantitative detail statis- 
tical moments for mass and for substructure samples whose 
densities vary by a factor of one hundred. 

On scales greater than of order a few Mpc, theoret- 
ical predictions are well represented in the point cluster 
version of the halo model. The point cluster model results 
range from a biased realization of halo correlations on large 
scales to intermediate scales, for which iV£ < 1, where sin- 
gle halo contributions dominate, but haloes are still unre- 
solved. In this limit the results of the halo model are in- 
dependent of details such as halo profile, asphericity, and 
concentration parameter. Intermediate results are robust; 
the variance £2 = fe/Nji steepens, approaching r~ 3 , an ef- 
fect seen in scaling studies (Hamilton et al. 1991; Peacock 
& Dodds 1996); and the amplitudes Sk are constant, as in 
the plateau seen in scale free models by Colombi, Bouchet, 
& Hernquist (1996). These results are independent of pro- 
file shape or bias. The halo model with Press-Schechter or 
Sheth-Tormen mass function allows us to compute from first 
principles values for the hierarchical amplitudes Sk for scale- 
invariant models with power spectrum P ~ k n . As shown 
in Figure 7, the halo model predictions are sensitive to the 
mass function and to mass cutoffs. For scale- free models with 
initial spectrum P(k) ~ k n , the halo model reproduces the 
general trends of Sk(n). Disagreements for more negative 
n are probably an indication of the difficulty of simulating 
these spectra. 

The largest halo has r s ~ 500 h~ kpc; on scales smaller 
than this, we must include the effects of finite halo size. Re- 



solved haloes introduce scale-dependent form factors in the 
fik,i(R), as in Section 4. Analysis of resolved haloes is made 
substantially more efficient by analytic expressions for the 
form factors, contained in Appendix C. Small-scale results 
suggest that the profile shape is different for mass and for 
substructures. In our limited efforts we have not found a 
profile shape that allows us to fit the shape of Sk on all 
scales in the resolved halo regime. For that matter, we do 
not really know that a universal profile shape applies for 
the distribution of galaxies within haloes of different size; 
A possible explanation is that tidal disruption leads to no 
universal profile that applies on all scales; or, the halo model 
picture itself may be oversimplistic. Nevertheless, on large 
scales our simulation and model results seem to have the 
potential to agree with observations (Ross et al. 2006). 

As observations become more and more precise, so it is 
increasingly important to be able to model clustering statis- 
tics with precision. This appears to be possible for mass on 
both large and small scales. On large scales, perturbation 
theory (biased linear theory) is accurate to better than 1%. 
On small scales, where statistics are dominated by tightly 
bound, high-density collapsed haloes, using published forms 
for the mass function f(m) and the concentration parame- 
ter c(m) with no attempt to optimize, the halo model repro- 
duces the variance for our simulation again to within a few 
percent. This is suggestive but not in itself a proof of the 
halo model; history has shown that there may be many con- 
structs that lead to the same two-point function; thus, it is 
a nontrivial result that the halo model also reproduces with 
accuracy the higher order correlation functions on small and 
large scales as well. There are somewhat larger deviations on 
intermediate scales, where the halo model predictions are 
too large by 5-10%, in a direction that is only made worse 
by including higher order perturbative corrections. This is 
a regime where the halo model seems to be least likely to 
be valid, where there is a significant amount of inhomoge- 
neously clustered mass not contained in spherical haloes; 
another interesting possibility in this regime is renormalized 
perturbation theory (Crocce & Scoccimarro 2006). The halo 
model predictions also match very well number count statis- 
tics on both large and intermediate scales, the point cluster 
regime, across all the different subcatalogs with different 
mass thresholds. This is perhaps not a surprise, since there 
is no "background" population of objects outside haloes; ob- 
jects in haloes account for all the objects there are. However, 
such precision does not seem to be within reach on small 
scales. The results we present use an isothermal profile with 
fixed c = 10, but this is at best only a first approximation. 
With a small number of haloes, we do not know the pro- 
file shape, although it seems that the NFW profile does not 
work, and we do not know how the halo radius or concen- 
tration parameter depends on mass. This may be the result 
of using substructures instead of galaxies in a full hydro- 
dynamic simulation; substructures in high density regions, 
may be tidally disrupted (Weinberg, Colombi, Dave & Katz 
2008). 

Halo model statistics computed over mass and number 
distributions taken from the simulation work well. It is in 
principle possible to compute correlations from first princi- 
ples, starting with a primordial power spectrum, using the 
Sheth-Tormen halo mass function and a prescription such as 
a Poisson satellite number. Application to scale-free Simula- 
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tions with initial spectrum P ~ k n gives plausible results for 
Skin), at least for n not too negative, once taking into ac- 
count finite simulation resolution. In practice, the predicted 
relative bias factors b g /bh do not quite match the numerical 
results, but this is probably due to finite volume effects. In 
particular, the halo five-point function is barely detected. 

In the end, on small scales there are substantial differ- 
ences between the discrete statistics of number counts and 
the continuous statistics of of mass. The distribution func- 
tion of halo occupation number has a behavior different from 
that of the distribution of halo mass, and factorial moments 
of discrete counts behave differently than moments of mass, 
even if the mean occupation number itself were a faithful 
tracer of total mass, all which contribute to differences in 
Sk, both in value and in shape as a function of scale, to the 
extent that it is not clear that the concept of bias between 
galaxy and mass statistics, even a non-linear bias, is a useful 
concept. 

It may sometimes seem that with a halo profile shape, 
mass function, concentration parameter, and asphericity all 
to be specified, the halo model is infinitely adjustable. How- 
ever, on intermediate and large scales much of this freedom 
disappears, and the model depends only on the compound- 
ing of statistics. In the halo model calculation, we see that 
the overall size of the correlation function £& or the ampli- 
tude Sk is determined by moments p, k — ( m k \ of mass or 
factorial moments ]lk = \ JV ) of halo occupation number; 
while details of shape on small scales provide information 
on the halo profile, ( F k \ . That the model can reproduce in 
detail the measured Sk for k — 3-5, simultaneously for both 
mass and number, and can handle probabilities as well as 
moments, is a nontrivial success. 
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APPENDIX A: EMPTY HALOES 

We show in this Appendix that for every composite distri- 
bution Pn that includes empty haloes with a nonzero proba- 
bility go 7^ 0, there is another with q' — that produces the 
same Pn- Thus, excluding (or including) empty haloes does 
not impose a physical restriction on the resulting occupation 
distribution Pn- 

Suppose we start with a distribution with qo ^ 0. Then, 
the total count probabilities will include contributions from 
many clusters with no occupancy, 

Po =Po+ Piqo +P'2q k + P:iq'o H (Al) 

-Pi = piqi + 2p 2 qiqo + 3p3<7i<?o H (A2) 

P2 =Pi<72 +P2(qi + 2q 2 qo) + p-i(3q 2 qo + Sqfqo) H (A3) 

Pi = Piq-i + P2(2qiq-2 + 2g 3 go) 

+P3(q( + QqiQ2qa + 3g 3 go) H (A4) 

We can easily create an occupancy distribution with no 
empty haloes while maintaining the same relative probabil- 
ities by defining a new set of probabilities q' n such q = 
and q' n — <?n/(l — qo) for n > 0. This distribution has the 
generating function 



g'i( z ) 



9i(z) -g 
1 - go 



(A5) 



Note that since <?i(0) = go, this gives g^(0) = 0, and §[{!) = 

3,(1) = 1- 

We can then modify the halo occupation number prob- 
ability p' n in what turns out to be a sensible way to produce 
in the end the same Pn- With q =0we must have Po = p'o, 
so we take 



Po = y^Pngp = gh(qo)- 

71 = 

Next, to have Pi = Piq[ 

. 00 

Pi9i 



Pi9i 



1 - go 
we take p[ = (1 



En-l 
np„gig 

n=0 



91 



dgh(z) 



Az 



(A6) 



(A7) 



qo)(dgh/dz)\ qo . Similarly, to have 



P2 =Piq2 +P2Qi 
= <?2 ^ np„q - 



En-l , 2 1 V^ I i\ ii-l 



+ gig } n{n- l)p„qo~ 



(A8) 



we take p' 2 — |(1 — go) 2 (d 2 (;j,/d2: 2 )| qQ . These p' n follow from 
the generating function 



g'h(z) =9h[(l -qo)z + qo\- 



(A9) 



The coefficient of each term in the expansion is a sum of 
products of positive numbers with positive coefficients, and 
so p' n ^ 0; and g' h (l) — gh(l — qo+c/o) — 1, so that each term 



must satisfy p' n ^ 1 and the distribution is properly normal- 
ized. The composition of these two modified distributions 
then gives 

G'(z) = g' h [g'i(z)] = g h [(l - qo)g'i(z) + qo] = gh[gi(z)], (A10) 

and so the same Pn, as desired. 

For the case of a Poisson cluster number distribution, 
Finkelstein (1988) shows that the revised p' n again arc 
a Poisson distribution with mean N' = A(l — go). The 
general case has essentially the same interpretation. The 
continuum (discreteness corrected) moments, generated by 
M(t) = G(l + t), follow from 

M' h {t) = G'(l +t) = G[(l - g„)(l +t) + go)} 

= G[l + (1 - qo)t] = M[(l - qo)t]. (All) 

Thus, continuum moments are scaled by a factor (1 — go)" 
which absorbed in the mean N' t = Ni{\ — qo), leave the 
correlations Jx n unchanged: the q' k are a discrete realization 
of the same underlying number density field. In a sense, this 
is the equivalent of including an unclustered background, as 
in equation (36). 

In general go can be mapped to any value < a < 1 by 
the transformations 



Ji(z) 



9i(z) 



g'h{z) =flh[(l -a)z + a}. 



(A12) 



and it remains true that the generating function of to- 
tal count probabilities is unchanged, G(z) — g'h[g'i(z)] = 
9h[gi{z)]- 



APPENDIX B: ALGORITHMS FOR 
COMPUTING THE COUNT-IN-CELL 
DISTRIBUTION FUNCTION 

In this Appendix we detail how the count-in-cell distribu- 
tion function Pn{(.) is estimated in these samples. There 
exist many efficient ways to measure this function in cubi- 
cal cells (e.g., Szapudi 1998; Szapudi, Quinn, Stadel &: Lake 
1999; Blaizot et al. 2006). The problem is however more 
intricate for spherical cells of radius £, which we prefer to 
use in this paper, since the analytical calculations are much 
easier to derive for these latter. Although it is rather usual 
and fair to approximate spherical cells with cubical cells of 
same volume with a small form factor correction (e.g., Sza- 
pudi 1998, and references therein), we prefer here to avoid 
this approximation. Then, the two most common ways of 
measuring function Pjv(£) for spherical cells are 

(i) The FFT method: it consists in assigning the particles 
to a grid of size -/V gr id using e.g. nearest grid point or cloud- 
in-cell interpolation (e.g., Hockney & Eastwood 1981), Fast 
Fourier Transform (FFT) the corresponding density distri- 
bution, multiply the result by the Fourier transform of the 
top hat filter in Fourier space, and then Fourier transform 
back to estimate function Pn(£)- Obviously the FFT method 
is valid only if the cell size is much larger than the size of a 
mesh element. 

(ii) The direct assignment method: it consists in creating a 
list of candidate cells positioned on a regular pattern of size 
iVg r id, then on scanning the list of particles and assigning 
them to each cell when relevant to augment the correspond- 
ing count. This method does not suffer the defects of the 
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Table Al. Parameters used to perform the count-in-cells measurements. The first line gives the inverse scale in units of the simulation 
box size L\, Q ^/£; the smallest and the largest scales, I/box/^ = 1024 and L\, a ^/£ = 8, correspond to £ = 0.2 h~ Mpc and £ = 25h~ 1 Mpc, 
respectively. The second line gives the size of the grid of sampling cells, N SI ^, used to perform the measurements at a given £ for the 
full dark matter sample, RAMSES; -/V grid = 2048 means that 2048 3 cells were used, corresponding to a minimum possible value of Pjv °f 
the order of 1.16 X 10 — 10 . The third line identifies the count-in-cell measurement method used for each scale under consideration, T for 
oct-tree walk, F for FFT, and D direct assignment. The fourth and fifth lines give -/V grid and the method for all the other samples. 
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FFT, so can be used even for very small scales, but can be- 
come prohibitive at large scales, when the network of cells 
starts to significantly self-overlap (i.e. a particle is assigned 
to a large number of cells) . 

Naturally, a good choice consists in using the direct assign- 
ment method at small scales and the FFT method at large 
scales. However, the very large number of particles, 512 3 , 
in our full RAMSES dark matter called for an additional al- 
gorithmic improvement, valid at intermediate scales. Our 
implementation, CountKD, is a code based on a decomposi- 
tion of space using standard oct-tree technique, similarly as 
in Szapudi, Quinn, Stadel & Lake (1999), but for spherical 
cells. The oct-tree decomposition allows one to approximate 
a spherical cell with a coverage of cubes of varying size, this 
latter decreasing (by factors of 2) when approaching the cell 
boundary. At some point, when maximum allowed level of 
refinement is reached (or when there is only one particle per 
oct-tree cell), the position of particles themselves is used to 
decide if they belong to the cell in consideration or not. Ob- 
viously, this method is efficient only if there are sufficiently 
many particles per cell, otherwise it is faster to perform di- 
rect assignment as explained above. For this work, we did 
not bother to find the optimal compromise between direct 
assignment, oct-tree walk and FFT method and performed 
the measurements as described in Table Al. The main goal 
was simply to reach a reasonable level of accuracy to sample 
correctly the large N tails of function Pn(€) in a reasonable 
amount of CPU time while avoiding as much as possible the 
FFT method which is sensitive to the pixelization of the 
data. 



APPENDIX C: FORM FACTOR INTEGRALS 

In this Appendix we present analytic expressions for the 
form factors F(x;y,c), the fraction of a halo with concen- 
tration parameter c and centre at r/r s = x that is contained 
within a spherical volume of radius R/r 3 = y. First, note 
that the volume V of the intersection of two spheres of radii 
a and b whose centres are separated by distance d is 

I(a, b,d), \a — b\ ^ d ^ a + b, 

0, d>a + b, 



V(a,b,d) = 



(CI) 



-^min(a 3 ,b 3 ), d ^ |a - b\, 



where the shared volume I (a, b, d) of spheres that partially 
overlap is 



I(a,b,d) = - 
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Next, note that a decreasing profile p(r) can be modelled as 
a sum of step functions, truncated at an outer radius r max . 
This means that the profile is composed of an ensemble of 
spheres of decreasing densities pi and increasing radii r,, 
i — 0,...,N, with r = 0, r N = r max . Then, for p(r) — p t 
for n-i < r ^ n, we can write 



p{r) 
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(C3) 



Now, let us compute the mass contained within radius R of 
the origin for a halo centred at a distance d, 



M(d,R) = / d A rp{d + r). 

•r<R 



(C4) 



From equation (C3), the calculation reduces to the sum of 
intersections of spheres with appropriate weights, 

M(d, R) = p(r m ^) V(R, r max , d) - ^ V(n,R, d)A Pi . (C5) 



In the continuous limit, we obtain the final expression 
M(d, R) = p(r max )V(i?, r max , d) 

I(r,R,d)^dr 
dr 
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(C6) 



We apply this in particular to the NFW and isothermal 
sphere profiles. 



CI Convolution of the NFW profile 

The truncated NFW profile writes, in scaled units 

PNFw(r) = J r(l + r)2' T ** °' (C7) 



r{l + r) 2 '' 
.0, 



r ^ c, 

r > c 
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The calculation of the integral gives 

M NFW (x) = - ( I(r, y, x)p(r) - ^£±12 ln (i + r ) 

x 



irr ivr 
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y 2 -(x + l) 

l + r 



r _ r _ 



- \ ^r 3 p(r) - 4tt — — + ln(r + 1) \ 
I A L r + 1 J J r=0 

-{|^ 3 p(0} r +p(c)V(c,2/,x), (C8) 

where p = pnfw, x = d/r s , Y" = R/r s , c = r ma x/r 3 , and 
r + = min(c, a- + y), (C9) 

r_ = min(c, |a; — j/|), (CIO) 

n = min[c, max(y — a;), 0)]. (Cll) 

Evaluated for y > c and a: = this gives the total mass, 

M = 47r [(l + c)ln(l + c)-c] (C12) 

(1+c) 

and the normalized NFW form factor is then 

47r[(l + c) ln(l + CJ - c]/(l + c) 

C2 Convolution of the Isothermal profile 

The isothermal profile writes, in scaled units 

r < c, 



Piso(r) = < l + r 2 ' " ' (C14) 

[0, r > c . 

The calculation of the integral gives 



— 2n(r — arctanr) 



Miso(x) = -h(r,x,d)p(r)-2n(r 
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-fWpwT +p(c)V(c,»,a;), (C15) 

where p = piso and r+, r_, and n are as in equations (C9)- 
(Cll). Evaluated for y > c and a; = this gives the total 
mass, 

M = 47r(c-arctanc), (C16) 

and the normalized form factor for the isothermal profile is 

F f~v^ Miso(a;;y,c) 

F IS o(a,Y,c) = Mc _ arctanc) . (C17) 



